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ABSTRACT 

We study the coorbital flow for embedded, low mass planets. We provide a simple 
semi-analytic model for the corotation region, which is subsequently compared to 
high resolution numerical simulations. The model is used to derive an expression for 
the half- width of the horseshoe region, x s , which in the limit of zero softening is given 
by x s /t p — 1.68(g//i) 1 / 2 , where q is the planet to central star mass ratio, h is the disc 
aspect ratio and r p the orbital radius. This is in very good agreement with the same 
quantity measured from simulations. This result is used to show that horseshoe drag 
is about an order of magnitude larger than the linear corotation torque in the zero 
softening limit. Thus the horseshoe drag, the sign of which depends on the gradient of 
specific vorticity, is important for estimates of the total torque acting on the planet. 
We further show that phenomena, such as the Lindblad wakes, with a radial separation 
from corotation of ~ a pressure scale height H can affect x S) even though for low-mass 
planets x s <C H . The effect is to distort streamlines and to reduce x s through the 
action of a back pressure. This effect is reduced for smaller gravitational softening 
parameters and planets of higher mass, for which x s becomes comparable to H. 

Key words: planetary systems: formation - planets and satellites: formation. 



1 INTRODUCTION 

Immediately after the d iscovery of the first extrasolar planet 
l|Mavor fc Quelozl Il995l ) , a Jupiter-mass planet in a very 
close orbit, it was realised that this class of planets, the 
Hot Jupiters, could not have been formed at their present 
location. In stead, they should have formed further out in 
the protopl anetary disc, and migra t ed inward a fterwards. As 
outlined in lGoldreich fc Tremainel l| 19791 . Il98(il ). planets em- 
bedded in protoplanetary discs indeed will undergo orbital 
evolution through disc tides, and a great deal of theoret- 
ical work has been dedicated to understand the direction 
and magnitude of plan etary migration (for an overview see 
iPapaloizou et alj|2007h . 

One can distinguish three types of migration. High mass 
planets, comparable to Jupiter, open up deep gaps in their 
discs, after which they migrate on approximately a viscous 
time scale dLin fc Papaloizou! Il986l ). This is called Type II 
migration l|Wardl 1997 ). Less massive planets, comparable to 
Saturn, embed ded in massive discs may undergo fast Type 
HI m igration l|Masset fc Papaloizoul 120031 ; iPepliriski et al.l 
2008a). Both Type II and Type III migration may be di- 
rected inward or outward, depending on local conditions in 
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the disc i|Crida fc Morbidellil [20071 : IPepliriski et al.ll2008bh . 
but the general trend is inward migration. 

Low-mass planets, up to a few tim es the mass of the 
Earth (M e ), undergo Type I migration l|Wardlll997T ). This 
type of migration is driven by a linear wave response in the 
disc, leading to a chara cteristic two-armed spiral pattern 
jQgilvie fc Lubowll2002h . The waves can be understood to 
be ex cited at Lindblad resonances (|Goldreich fc Tremainel 
Il979l ) and lead to a torque on the planet that is due 
to asymmetries in the density , pressure and rotation pro- 
file in the disc l|Wardl I1997T ). The resulting migration 
direction is inward for all reasonable disc p arameters 
jKorvcanskv fc Pollacklll993l ; iTanaka et alj|2002h . 

Apart from this wave, or Lindblad, torque, em- 
bedded planets are also s ubject to corotation torques 
jGoldreich fc Tremaine1ll979h . Two descriptions of the coro- 
tation torque exist in the literature. As advocated by 
iGoldreich fc Tremainel (| 19791 ). one can perform a linear anal- 
ysis of the corotation resonance, which leads to a torque 
proportional to the radial gradient of specific vorticity, or 
vortensity, in the unperturbed disc. Semi-analytical and nu- 
merical studies lead to expressions f or the total corotat ion 
torque (|Korvcanskv fc Pollacklll993l ; ITanaka et al.ll2002f ) in 
two- and in three-dimensional discs. We will refer to this 
torque as the linear corotation torque. 

A different view on the corotation torque was given by 
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IWardl(|l99ll ). who considered the torque due to material near 
the orbit of the planet that executes horseshoe turns. The 
total corotation torque is again proportional to the radial 
gradient of specific vorticity in the unperturbed disc, but 
the model contains a free parameter x s , the width of the 
horseshoe region. We will refer to this torque as the horse- 
shoe drag. The relation between the two descriptions has 
never been clarified so far. 

Linear theory has been compared successfully 
against numerical h y drodynamical simulations, in 2D 
llD'Angelo et al.l |2002|; iNelson fc Papaloizou! [20041 ) as well 
as in 3D l)D'Angelo et al.l l2003h . At this point we remark, 
however, that in these studies, only discs with small 
gradients in specific vorticity were considered, that is, cases 
where the corotation is supposed to be weak. For shallow 
surface density gradients, resulting in strong corotation 
torques, intermediate mass planets m ay reverse their 
direction of migration ijMasset et al.ll2006h . 

All studies mentioned so far have made use of the sim- 
plifying assumption that a barotropic (or isothermal) equa- 
tion of state applies, in which case no energy equation needs 
to be solved. The dram atic effects of releasing this assump- 
tion were first noted in iPaardekooper fc Mellemal (|2006al l . 
where it was shown that low-mass planets can migrate out- 
ward in non-isothermal discs. Subsequently it was realised 
that this was d ue to a radial entropy gradien t in the un- 
perturbed disc (jPaardck ooper fc M cllcma 2008), giving rise 
to a s trong, positive corotation torque. iBaruteau fc Masse! 
(2008) provided a linear analysis of the problem, and ar- 
gued that a linear effect due to a background entropy gra- 
dient can be strong en ough to reverse the torque on low - 
mass planets. However, IPaardekooper fc Papaloizou! (|2008l ) 
showed that this linear effect is in fact small, and that a non- 
linear effect is responsible for the torque reversal as seen in 
IPaardekooper fc Mellemal |2006al ). 

The non-linear t or que, as studied in 
IPaardekooper fc Papaloizou! l|2008l '). is closely related 
to the idea of horseshoe drag originally introduced by by 
IWardl ([ 199ll ) as both of these are produced by disc material 
undergoing horseshoe turns in the neighbourhood of the 
planet. For this reason it is very important to have a clear 
understanding of the horseshoe drag and in particular its 
relationship to linear corotation torques that are often 
used to estimate torques arising from coorbital effects. In 
this paper, we will provide an analysis of the horseshoe 
region for low-mass planets, in particular we determine its 
half-width x s for softening lengths ranging between zero 
and the disc scale height. For simplicity we shall return 
to considering a barotropic equation of state and work 
with a two dimensional model. Since the horseshoe drag 
is proportional to x$ (see IWardlll99ll ). a good estimate of 
this parameter is critical in determining the total torque 
on the planet which is found in to be significantly larger 
than estimates based on linear theory for zero softening. 
In an accompanying paper, we perform a general study on 
the behaviour of the torques and their dependence on other 
parameters such as the disc viscosity. 

The plan of this paper is as follows. In Section [2] we 
review the basic equations and introduce our local model. 
In Section^ we study the structure of the corotation region, 
and in Section|4]we analyse the resulting streamlines. These 
are then compared with numerical hydrodynamic simula- 



tions in Section [5] after which we give a brief discussion and 
conclusions in Section [6] 



2 BASIC EQUATIONS AND DISC MODELS 

The basic equations are those of the conservation of mass, 
momentum and energy for a two dimensional disc in a frame 
rotating with angular velocity Q p . we adopt a cylindrical 
polar coordinate system (r, ip, z) with origin (r = 0) located 
at the central mass. The disc then occupies the plane z — 0. 
The continuity equation and the equation of motion take 
the form 

rlY 

^ = ~V.(Ev) (1) 
and 

Dv - 1 

— + 2fi p kx v = --Vn-V<E> (2) 

respectively. Here, E denotes the surface density, v = 
(v r , v v ) the velocity, k denotes the unit vector in the vertical 
direction, and n the vertically integrated pressure. Thus 



n 



Pdz. 



(3) 



The total potential $ is taken to be $ = <E>G + ^pV 2 /2, where 
<E>G is the gravitational potential. The convective derivative 
is defined by 

We adopt a barotropic equation of state such that 

H = F(E), (5) 

with F(S) being a prescribed function of E. The square of 
the sound speed is given by 

2 _ dF(E) 
dE ' 

When a power law is adopted such that 
F = A"E", 



(6) 



(7) 



with K and /? being constants, c s = /3(n/E). The discs 
we consider are assumed to be low mass with the Toomre 
parameter Q — (Qc a )/(nG'E) ^> 1. The self-gravity of the 
disc is accordingly neglected. 

2.1 Global studies 

The gravitational potential is assumed to be due to the cen- 
tral mass and perturbing planet such that $g = $go + $Gp, 
when the z dependence is neglected. These are given by 

-GM, 



$go = 
and 
$G P = 



(8) 



—GM P 



+ 



yjr 2 + 7-2 — 2rr p cos(Ai^) + b' 2 r 

GM p rcos(A.ip) 



(9) 



where Atp = ip — <p p . In the above M» denotes the mass of 
the central object, with M p , r p , and tp p denoting the mass, 
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orbital radius and angular coordinate of the protoplanet re- 
spectively. The gravitational softening parameter is b. The 
last term in (|9} is the indirect term which accounts for the 
gravitational acceleration of the origin of coordinate system 
due to the action of the protoplanet. 

Rather than neglecting the vertical dependence of the 
perturbing potential, one may adopt a vertical average based 
on the density, p, thus one replaces the expression ([9} by its 
vertical average 



p^Gpdz. 



(10) 



This procedure is more consistent with the view that the 
two dimensional disc representation should be derived from 
applying a vertical averaging procedure and accordingly pro- 
vides a vertical average. On the other hand if vertical mo- 
tions near the midplane are small, use of the unaveraged 
form may be appropriate for rep resenting the struct ure there 
as in a stacked layer model (eg. iMasset et al.1120061 ). 

When Q p — the reference frame is non rotating 
but non inertial as the origin accelerates together with the 
central mass due to the action of the protoplanet. Nu- 
merical calculations are most conveniently performed in a 
frame corotating with the protoplanet. Then Q p becomes 
the circular Keplerian angular velocity at radius r p . At a 
general radius, r, the Keplerian angular velocity is given 
by n K (r) = (GM«/r 3 ) 1/2 . Thus D. K (r p ) = Q p . The lo- 
cal vertical scale height of the disc is then defined through 
H — c s /Qk- Hydrostatic equilibrium in the z direction im- 
plies that 



po exp 



(z'/H 2 )dz 



(11) 



where po is the midplane density. The value of H at the disc 
midplane is a measure of the local disc semi-thickness. From 
now on H will be stand for this quantity. 



2.2 Local Description 

As we are interested in a local region close to the planet, 
it is useful to work with a form of the basic equations 
adapted for thi s purpose. W e use the well known shearing 
box formalism (|Wardlll99ll ). This uses a Cartesian coordi- 
nate system (x,y,z) corotating with the protoplanet and 
with origin at the centre of the planet. The disc velocity in 
the unperturbed state with no protoplanet, is in the local 
approximation^ = (v x ,v y ,0) = Vo = (0, — 3f2 p :r/2, 0), cor- 
responding to a linear shear. When a protoplanet is present 
the equation of motion may be written 



9v 



+ v ■ Vv + 2Sl p k x v 



:Vn-v$ L , 



(12) 



where <J>l = 5>Gp — 3£l p x 2 /2. The protoplanet potential is 
here taken to be given by 



$G P = 



-GM P 



sjx 2 + y 2 + b 2 rp 



(13) 



The indirect term being small compared to the direct poten- 
tial close to the planet is neglected. The continuity equation 
remains of the same form as equation JT}. 



3 THE COROTATION REGION 

We now consider the structure of the corotation region. We 
begin by considering steady state solutions of the basic equa- 
tions. In this case the two components of the equation (|12|l 
are 



v • X7v x — 2Q p v y — — 



Ox 



and 



v • Vv y + 2Q. p v x 



dy 



(14) 



(15) 



where we have introduced the enthalpy, tu(£), which is de- 
fined through dw/dT, = cf/E together with the specification 
that w(0) = 0. 



3.1 A simple one dimensional model 

We now simplify the problem by in the first instance ne- 
glecting the degree of freedom corresponding to epicyclic 
motions. These are not expected to play a major role in the 
horseshoe region. One expects a balance between the Corio- 
lis force and potential gradient with v • Vv x being negligible 
in equation (|14[) . Neglecting this term is precisely the ap- 
proximation often used in celestial mechanics to enable the 
derivation of a second order differential equation describing 
particle motion on horseshoe orbits. It breaks down only 
close to the protoplanet. In the next subsection, we discuss 
a more complete description of the coorbital region. 

Adopting this approximation and using equation (|14[l 
we set 



1 

2U 



dx 



(16) 



where \ — w + with ip being the stream function. Equa- 
tions (|14[) and (|15l) may also be written in the form 



2fi p + 



&Vy 

dx 



= -V 



1 2 

-Vy + W + "PL 

2 " 



(17) 



From this and the steady state form of the continuity equa- 
tion fl]) it follows that both 



E oW = [ -jVy + W + $L 



and £o(ip) 



2fi D + 



dl ' ll 

dx 



are constant on streamlines. These are statements of the 
conservation of the Bernoulli constant and the specific vor- 
ticity or vortensity on streamlines under the approximation 
used here which has the consequence that the contribution 
from the radial velocity is neglected. We note that the func- 
tional forms of Eo(tp) and £o(V>) cannot be determined from 
the inviscid equations but have to be prescribed externally. 
From (1161) it then follows that 



1 



dx 2 + p 



(18) 



Using the fact that tu(£) = x ~ to eliminate E, a single 
equation for \ then results. We note that when a power law 
equation of state is used with /3 = 2, see iff}. «;(£) = 2KY*, 
and we have the very simple relation 



(19) 
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We thus obtain an equation for \ °f the form 



1 d\ 



(20) 



(21) 



2fi p dx 2 

Setting X = Y - 3fi 2 :r 2 /2, this leads to 

y = cb | j Q2Y i ^ P 
Gp 2CoSfi p da; 2 2£ E' 

We comment that although we adopted a power law equa- 
tion of state with [3 — 2 to obtain (|21[) . because w is a linear 
function of E in this case, the same equation applies to any 
barotropic equation of state in the regime where the surface 
density variations are small enough that only first order vari- 
ations in the equation of state need to be considered. Then 
c 2 /E is taken to be the constant background value. Of course 
when (3 = 2, c 2 /£ = 2K = constant regardless. 

3.1.1 Solution for the coorbital region 



Equation (|21[) contains no y derivatives and can accordingly 
be solved as a second order differential equation. Adopting 
the boundary condition that Y is bounded for \x\ — > oo the 
solution can be written down in the standard form 



Y = 



|°° G{x,x) (n 2 p + * Gp (x',y) 



dx\ 



(22) 

where unless explicitly stated, quantities in this and other 
integrands are evaluated at (x' ,y') and the Green's function 
satisfies 

d 2 G(x,x') 



G{x,x ) 
ox' 2 



/2|o|np\ =s{x _ x y (23) 

V s / x=x' 

with G(x,x') — > for \x'\ — > oo. Having found Y = x + 
3f2pX 2 /2, the streamlines can be found from 

which is constant on them. 



3.1.2 Constant vortensity 

The corotation region is generally of small radial extent so 
that the quantity £o> representing the variation of vorten- 
sity should be approximately constant unless the profile is 
very sharp. Therefore for reasonable smoothly varying cases 
vortensity variation is not expected to have much effect on 
the the coorbital structure except possibly for large soft- 
ening cases perturbed by the density wakes associated with 
Lindblad torques (see below). Accordingly we shall specialise 
to the case of constant vortensity while retaining /3 — 2, so 
that both E/Cg and £o are constant. By considering large dis- 
tances from the protoplanet, we infer that £o = fi p /(2Eo), 
where Eo is the uniform surface density at large distances 
from the protoplanet. Determination of the Green's function 
is straightforward in this case. We obtain 

G(x,x') = -±-exp(-\k\\x-x'\), (25) 



where 

O 2 

k 2 = )}v= H - 



where c s o is the uniform sound speed at large distances from 
the protoplanet and H is the scale height. 

Although solution of (|21[) is straightforward, we com- 
ment that the scale of the decay of the Green's function 
is the scale height and thus phenomena that distance away 
from the protoplanet and not included here, such as the 
prominent wakes associated with Lindblad torques, can dis- 
tort the flow (see below). In principle this effect could be 
included through the boundary conditions on G but we shall 
not consider it further in this section. 



3.1.3 Asymptotic series and softened gravity model 

We note also that in this limit, provided £o also varies on a 
length scale significantly exceeding, H, one can find a solu- 
tion of equation (|21|) in the form of an asymptotic series in 
ascending powers of Cg or equivalently H 2 /\y\ 2 . The zero or- 

c 2 f2 

der solution being simply Y = "^Gp + 2 | s , and we note here 
that the second term is a constant if £o is constant and may 
be discarded, or equivalently \ — "^l- The streamlines are 
then found from (|24[) using this value of \. We comment that 
this result would be obtained if pressure was neglected com- 
pletely. Thus we describe it as the softened gravity limit and 
it should apply to streamlines at a distance greater than ~ H 
from the protoplanet. It can also be obtained when the ap- 
proximation of neglecting the acceleration in the x direction 
is not made at the outset. We now consider a modification 
of the above model that enables a more realistic treatment 
of regions close to the planet. However, this makes it two 
dimensional and accordingly more complex. 



3.2 Modification close to the planet 

We begin by noting that the system may be regarded as 
being governed by the components of the equation of motion 
(|14|) and (|15[) together with a complete statement of the 
conservation of vortensity in the form 

=(*>, + (27) 

This equation should be regarded as replacing the continu- 
ity equation in those governing the model. The simple one 
dimensional model is obtained by using only (114[) and (|27[) 
with v x being neglected. Equation (|15[) is then used as an 
auxiliary equation to subsequently determine v x . We now 
retain v x in equation (|27[1 which now needs to be specified 
using equation (|15[) in advance. We shall retain the approx- 
imation of neglecting v x or equivalently the radial accelera- 
tion equation in (|15p . Thus as before we have 



i d x 

Vy 2fi p dx ' 
Equation (|14|l gives 



r, ! 2<> ; , + ^j- 



0(x + t;g/2) 
dy 



(28) 



(29) 



(26) 



We make the approximation of replacing v y in this equation 
by the unperturbed Keplerian value v y = — 3xn p /2.This 
leads to 

V X = -^. 

fi p dy 



(30) 



On the corotation region for low-mass planets 5 



We comment that we have found the above equation to be 
satisfied close to the planet in numerical simulations (see 
below). 

Using (|28l ), (|30l) and again adopting a power law equation 
of state with /3 = 2 and using (|19[) . equation (|27|l leads to 



p2y= 2^|oE (F _ $Gp) 



where the operator 



d 

dx 2 



<9y 2 



(31) 



(32) 



and we recall Y = x + 3f2 2 :r 2 /2. The required solution of 
equation (|31[) can be written as 

Y = -f° G(x, y, x', y') (q 2 p + (^p) *g p ) dx'dy', 

(33) 

where G(x,y,x' ,y') is a two dimensional Green's function. 
When the vortensity is constant, and the Green's function 
that vanishes at oo, it is given by 



G(x,y,x',y') = - — K (ky/(x - x') 2 + (v - j/) 2 /4), (34) 

where as in standard notation Tfo is the Bessel function 
(this is readily obtained after scaling the y coordinate to 
transform TD 2 into the Laplacian). When comparing with 
the corresponding result obtained from the one dimensional 
treatment given by equation (122[) , we see that there is a two 
dimensional as compared to a one dimensional integration. 
The additional integration over y' acts like a smoothing on 
a length scale H. 

The streamlines are in general determined from the 
Bernoulli condition that (v 2 /2 + w + $l) be constant on 
them. This reduces to the condition given through equation 
(|24|) under the approximation scheme used to derive (|31[) . 
We further remark that the scale of the flow in the y di- 
rection then becomes larger than H so we may perform the 
integration over y keeping other quantities fixed. Using the 
result that 
1 
2 



K (ky/(x - x') 2 + (y - y 1 f r /i)dy' 



exp(— \k(x — x)\ 



(35) 



we then recover equation (1221) for the simple one dimensional 
model. 

We may also develop an asymptotic sequence as in sec- 
tion |3T3] starting with Y = $ Gp - fip + c 2 Sl p /(2£ E), with 
the last two terms being constant, and then iterating (|31fl to 
find successive corrections to Y. Thus use of the two dimen- 
sional model imparts a smoothing in the y direction, with a 
length scale H, to the one dimensional model. This causes 
differences close to the planet for small softening. 



4 STREAMLINE CALCULATIONS 

We have used the simple one dimensional model to obtain \ 
and determine streamlines for the constant vortensity case 
using a range of softening lengths. We perform the integral 
specified in equation (|22[) in conjunction with equation (124[) . 
For comparison purposes we also obtained streamlines using 



equation (|33[) . instead of (|22[) . In this case the Green's func- 
tion and integration are of course two dimensional. 

4.1 Dimensionless scalings 

Adopting H evaluated at the origin for the unperturbed 
flow as the unit of length, it is straightforward to see that 
apart from £o, the distribution of specific vorticity, if b/h 
is fixed, the streamlines are characterised by only one pa- 
rameter, q/h 3 , where q is ratio of the mass of the proto- 
planet to the mass of central star and h = H/r v is th e disc 
aspect ratio (see e.g. iKorvcanskv fc Papaloizoulll996l ). We 
work with the power law equation of state with /3 = 2. In 
that case c 2 /E is constant. In the local model, the adop- 
tion of a linear shear means that the background vortic- 
ity is a constant = (l/2)fi p . Therefore for constant vorten- 
sity the background surface density is a constant = Eo and 
£o = fi P /(2Eo). The units are arbitrary and so we choose 
these such that E = 1 at the origin in the unperturbed flow. 

Streamlines for constant vortensity obtained for the 
simple one dimensional model from ([22} with q/h 3 = 0.0252 
and b/h = 0.6 are shown in Fig. [T] This parameterization 
corresponds to 1 Mq in a disc with H/r v ~ 0.05. The 
streamlines fall naturally into two groups, the first coming 
from y > and x > and the second coming from y < 
and x < 0. For each of these classes a subset passes by 
the planet while the remainder undergoes a horseshoe turn. 
Those undergoing horseshoe turns constitute the horseshoe 
region, which is separated from the remaining domain by 
two separatrices. For the constant vortensity case, the flow 
has both left-right and up-down symmetry implying an X 
point at the centre of the protoplanet. When the vortensity 
is not constant this symmetry is in general lost. 

The streamlines shown in the upper left panel of Fig. 
[1] corresponding to q/h 3 = 0.0252 and b/h = 0.6. have a 
horseshoe width of ~ 0.25H. To illustrate the dependence on 
q we also show the streamlines for q/h 3 — 0.0126 and b/h = 
0.6. These suggest that the horseshoe width is proportional 
to q 1 / 2 . This result can also be derived by considering the 
streamline at the centre of the protoplanet. At this location, 
where x — dx/dx — 0, equation (|33[) gives 



X = Y 



1 



G(0, 0, x',y)$Gpdx'dy' + c 2 fi p /(2£ E), 

(36) 



there being a corresponding expression derived from equa- 
tion (|22l) . The horseshoe width, x s , is obtained by equating 
this expression to c 2 f2 p /(2£oE) — 3Q,pX 2 /8, being the value 
of Eo(tp) obtained from ((24} at large distances. Hence very 
generally, 



(37) 



To be more precise, performing the integral (I36|l for the two 
dimensional case we find that 



K {r)r 



r% 3irhJ J ^ r 2 (l + 3sin 2 8) + b 2 /h 2 
which can be simplified to 



dOdr, (38) 



16g 



Ko(r)r 



E 



r% 3nhJ ^/4 r 2 + b 2/ h 2 \V 4r 2 + b' 2 / h 



3r 2 



dr, (39) 
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Streamlines Streamlines 




Figure 1. Streamlines for constant vortensity obtained using the simple one dimensional model. The plots are for q/h 3 = 0.0252 
and b/h = 0.6 upper left panel, q/h 3 = 0.0126 and b/h = 0.6, upper right panel, q/h 3 = 0.0252 and b/h = 0.3, lower left panel and 
q/h 3 = 0.0252 and b/h = 0.0252, corresponding to the softening parameter being equal to the Bondi radius, lower right panel. For this 
and other similar figures the unit of length is the disc scale height, H, evaluated at the origin of the unperturbed flow. 



where E denotes the complete elliptic integral of the first 
kind. For b = 0, this gives 




which gives the horseshoe width in the limit of zero softening 
for the two dimensional model. We also note that in the limit 



6> h, 




(41) 



which is the same result that would be obtained by neglect- 
ing pressure effects altogether. 

To investigate the dependence on the softening param- 
eter, we plot streamlines for the cases q/h 3 = 0.0252, with 
b/h = 0.3 and also with b/h = 0.0252 in Fig. □ The second 
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Streamlines Streamlines 




Figure 2. Streamlines for constant vortensity found using the two dimensional Green's function. The plots are for q/h 3 = 0.0252 and 
b/h = 0.0252, corresponding to the softening parameter being equal to the Bondi radius, upper left panel, and q/h 3 = 0.0252 with 
b/h = 0.1, upper right panel. In the lower left panel, q/h 3 = 0.0252 with b/h = 0.3, and in the lower right panel q/h 3 = 0.0252 with 
b/h = 0.6. The horseshoe region is wider for smaller softening lengths and converges to a well defined structure for b — > 0. 



case corresponds to the softening parameter being equal to 
the Bondi radius Cg/(GAf p ) for 1 M9 and a disc aspect ra- 
tio of H/r p — 0.05. Interior to the Bondi radius the gravity 
of the planet dominates and we expect a hydrostatic struc- 
ture that does not participate in the horseshoe dynamics. 
We here remark that equation (|22[) obtained for the simple 
one dimensional model diverges logarithmically for 6^0, 
but (|33p is convergent. 

In general as expected from (|37[) the horseshoe width 



increases as 6 decreases and the streamlines become more 
compressed and nearly horizontal near y — indicating that 
the acceleration in the radial direction should not be ne- 
glected for small softening as has been done in the simple 
one dimensional model. In spite of this, the behaviour of the 
streamlines for y > H is very similar for fixed q and the 
different values of b in these cases. This is in line with the 
existence of the asymptotic solution discussed above. 

Streamlines for constant vortensity when the two di- 
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Figure 3. Streamlines, obtained from hydrodynamical simulations 
right panels: b/h = 0.025. Top panels: full potential without cut off, 



mensional Green's function is used are given in Fig. [2] These 
are for q/h 3 = 0.0252 with b/h = 0.0252, corresponding 
to the softening parameter being equal to the Bondi ra- 
dius, together with plots obtained for b/h — 0.1, b/h = 0.3 
and b/h = 0.6. Use of the two dimensional Green's func- 
tion smooths the potential and makes the horseshoe re- 
gion narrower compared to the simple one dimensional case. 
When b = 0.6h, the horseshoe width was 0.211/ being about 
20% smaller than that found using the simple one dimen- 
sional model. However, the horseshoe width converges in the 
limit of zero softening being about 30% wider than when 
b/h = 0.6. 

We remark that using a potential vertically averaged 



0.96 0.98 1 bo 1.02 1.04 




x 



0.96 0.98 1.00 1.02 1.04 




x 



for constant vortensity and q/h 3 = 0.0252. Left panels: b/h = 0.6, 
bottom panels: the cut off procedure was adopted. 



with weight p, equivalent to a vertical smoothing, but re- 
taining the formalism leading to equation ([22]) has a similar 
effect to using the two dimensional Green's function includ- 
ing convergence for b — > 0. 

The discussion presented above was for low-mass plan- 
ets, for which the Bondi radius is smaller than the radius 
of its Hill sphere r p (q/3) 1 ^ 3 . This condition is equivalent to 
the requirement that q < /i 3 /\/3- In this regime we found no 
circulating streamlines close to the planet when the vorten- 
sity is constant. For more massive planets with q > h 3 /\/3, 
a circulating region is expected to interior to the H ill sphere. 
This l imit was recently considered analytically in iPepliriskil 
(2008). We show an example of the resulting streamlines ob- 
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x(R H ) 

Figure 4. Streamlines close to the planet in the limit q 3> h 3 /^/3, 
in which case a circulating region within the Hill sphere arises. 
The unit of distance in this figure is Hill radius Rjj = (g/3) 1 ' 3 r p . 
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Figure 5. Streamlines close to the planet for q/h 3 = 0.0252 
and b/h = 0.6 for a disc with constant vortensity. The planet is 
indicated by the filled circle, and the dashed contour indicates 
the Bondi radius. The cut off procedure was adopted for this 
simulation. 



tained from our model in this limit in Fig. [4] Even though 
important non linear effects such as the onse t of gap forma- 
tion are not included (|Papaloizou et al.l l2007) there is indeed 
a circulating region within the Hill sphere in this case. 

However such closed streamlines might in general be 
expected significantly interior to the Bondi radius for matter 
bound to a low mass planet. We comment that this issue is 
sensitive to model details such as the specification of the 
vortensity profile. To illustrate this, consider the situation 
when this bound matter appears almost stationary in the 
rotating frame such that circulating streamlines would be 
expected. This matter would have to have significantly lower 
vortensity than its surroundings. To see this we note that 
for a staionary solution x — 0, and Y — 3Q.pX 2 /2. With this 
specification, equation (|3ip gives 



2Q p cj _ 3npV 



Gp, 



(42) 



setting the vortensity profile to be such that £o = 2fi p /E, 
we obtain 



<f> 



Gp. 



(43) 



which is the expected condition for hydrostatic equilibrium 
of the protoplanet under centrifugal and tidal forces for our 
baratropic equation of state. Note that to achieve this, the 
surface density should rapidly increase and hence the pre- 
scribed vortensity should correspondingly decrease towards 
the center of the protoplanet. On account of the attainment 
of a limiting form of the horseshoe region once the softening 
lenght b is smaller than the Bondi radius, we do not ex- 
pect the details of behaviour at interior radii to significantly 
affect the results. 
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Figure 6. Streamlines close to the planet for q/h 3 = 0.0252 
and b/h = 0.6 for a disc with constant vortensity. The planet is 
indicated by the filled circle, and the dashed contour indicates 
the Bondi radius. The full planet potential was adopted with no 
cut off. 



5 NUMERICAL SIMULATIONS 

In this section, we compare the streamline calculations of 
Section [4] with fully non linear hydrodynamical simulations. 
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5.1 Set up 

We u se the RODEO method l|Paardekooper fc Mellemal 
2006b) in two spatial dimensions, on a regular grid which 
when at its most extended, runs from r — 0.5r p to r = 1.8r p 
and which covers the whole 2tt in azimuth. Since we want 
to resolve the horseshoe region for even the smallest plan- 
ets we consider, a relatively high resolution is used. For the 
most extended grid this has 1024 cells in the radial and 4096 
cells in the azimuthal direction. Then the resolution at the 
location of the planet is approximately 0.0015r p in both di- 
rections. Tests have shown that this resolution is sufficient 
to capture the horseshoe dynamics for x s > 0.004r p . We al- 
ways ensure that we resolve the softening length br p by at 
least 3 grid cells, which means that for the smallest values 
of b we consider an even higher resolution is was adopted. 
We take the disc to be inviscid and isothermal with uniform 
specific vorticity such that £ oc r~ 3 ^ 2 for the unperturbed 
initial state. 

We consider two kinds of simulation, the first adopted 
the perturbing potential given by equation © and adopted 
the most extended domain. The second type of simulation 
considers only the coorbital region extending from r p — 2H/3 
to r p + 2H/3 and occupying the full 2tv in azimuth. Non re- 
flecting boundary conditions are applied at the radial bound- 
aries such that material is allowed to leave and enter freely. 
We refer to this second type of simulation as having em- 
ployed a cut off procedure. Results obtained with this type 
of simulation were checked with simulations employing the 
most extended domain but modifying the protoplanet per- 
turbing potential such that it is given by equation ((9} for 
r — r p \ < 2H/3, and zero otherwise. These gave very similar 
results. The cut-offs applied in these simulations exclude the 
bulk of the contributions from the Lindblad torques to the 
flow. Since either the planet potential vanishes in the region 
where they are normally generated or that region falls out- 
side the computational domain. As remarked in Section |3j 
features on the order of one scale height away, such as phe- 
nomena associated with Lindblad torques, may affect the 
coorbital region even when this is much narrower than H. 
Differences between the results obtained from the two kinds 
of simulation we performed are consistent with this suppo- 
sition. 



5.2 Streamline analysis 

We start by comparing the flow on a scale of H to that 
found using the model of Section [3] in Fig. [3] for two differ- 
ent values of the softening parameter b and the two types of 
simulation. To obtain q/h 3 = 0.0252 we set up 1 around 
a Solar mass star in a disc with h = 0.05. We compare the re- 
sults illustrated in this figure with the corresponding results 
obtained using the two dimensional Green's function shown 
in Fig. [21 For the case of large softening, b/h = 0.6, without 
a cut off, although the streamline pattern appears similar, 
the width of the horseshoe region is about 22% smaller than 
that depicted in Fig. [21 When the cut off procedure is ap- 
plied, the horseshoe region increases in width by about 50% 
(see below). This indicates that phenomena located in the 
region where the Lindblad torques are generated may play 
a significant role in shaping the coorbital region. 

For smaller softening such that b/h = 0.025, corre- 



sponding to the softening parameter being equal to the 
Bondi radius, the differences between the two types of simu- 
lation is less extreme. The measured increase in a; s obtained 
on applying the cut off procedure amounts to 20%. In this 
case there is very good agreement between the simulation 
without a cut off and the model of Section with the two 
dimensional Green's function illustrated in Fig. [2] with the 
values of x s differing by several percent. 

The dependence on softening can be understood as fol- 
lows: the strength of phenomena related to Lindblad torques 
is determined by the planet potential at approximately a dis- 
tance H from the planet, while the width of the horseshoe 
region depends on | "3?Gp | at the location of the planet. For a 
softening parameter comparable to h, the planet potential at 
these two locations is comparable. Therefore, it is relatively 
easy for Lindblad torque related phenomena to affect the 
coorbital region. For smaller softening, [<&Gp| at the location 
of the planet increases, while the value at a distance H re- 
mains largely unchanged. For b <C h, we expect the effect of 
phenomena originating a distance ~ H from the planet to 
be smaller, and therefore there should be better agreement 
between the numerical simulations and the model of Section 

El 

The effect of introducing the cut off is to produce larger 
corotational speeds directed towards the planet. These are 
eventually slowed down as the pressure gradient reduces, 
producing a stagnation point at, or very close to the planet's 
location. The faster moving material can originate further 
from the planet and so is associated with an increased horse- 
shoe width. 

The weaker corotational flow that occurs without the 
cut off is associated with significantly increased pressure in 
the region near <p = <p p . This back pressure is affected by 
conditions a distance ~ H from the planet and it can distort 
the streamline pattern close to the planet so that it becomes 
asymmetric as shown in Figs. [5] and [6] For the larger soft- 
ening cases, the stagnation point is displaced azimuthally 
a distance ~ br p from the location of the planet and the 
horseshoe width is reduced (see Fig. [3]). 

However, when the cut-off is applied, there is a sin- 
gle stagnation point at ip = ip p , slightly displaced from the 
radial location of the planet due to the radial pressure gra- 
dient of the unperturbed disc. The latter is a very minor 
effect that is absent from the local models on account of 
their strict symmetry. The local models always produce a 
single stagnation point at the location of the planet, a sit- 
uation that is essentially recovered for small softening in 
simulations without a cut off (see also Fig. [3j . 

We take a closer look at this back pressure effect in 
Figs. [Jj and [8] where we consider a slice through r — r s t ag , 
where r sta g is the radial location of the stagnation point. 
Because of the radial pressure gradient of the unperturbed 
disc, r st ag 7^ r p (see Figs. [5] and [6]). From equations (|28[1 
and (|30[) we know that the velocity is directly related to the 
gradient of %■ First of all, we can check if this approximation 
is valid. From the left panel of Fig. [7] we see that for both the 
full and the cut-off potential there is very good agreement 
between both v x and v y as measured from the simulation and 
the values expected from equations (128[) and (|30[1 . Therefore, 
their assumption in the two dimensional model of section [3] 
should be a good approximation. 

Note, however, the strong differences between the black 
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Figure 7. Velocities, enthalpy and planet potential at the radial location of the stagnation point for q/b? = 0.0252 and b/h = 0.6. Black 
curves indicate that the full planet potential was used, grey curves indicate results obtained with the cut-off procedure adopted. 




y 

Figure 8. Derivatives of the enthalpy and the potential at the 
orbital radius of the stagnation point for q/h? = 0.0252 and b/h = 
0.6. 

curves for the full potential and the grey curves for the cut- 
off potential. Velocities are reduced due the back pressure 
originating at distances ~ H from the planet when there is 
no cut off. The right panel of Fig. [7] clearly shows a higher 
peak in the enthalpy for the potential without a cut off. 
When added to the protoplanet potential the total is less in 
magnitude in that case and thus we can see from equations 
(|28p and (|30p that smaller inflow velocities will accordingly 
be produced. 

The azimuthal shift of the stagnation point can be ex- 
plained as follows. At a stagnation point, we must have that 



dw/dy = —d&h/dy (see equation JTSJ)), since all velocities 
must vanish. When the cut-off is applied, there is only one 
stagnation point possible (see Fig. [8} at the azimuth of the 
planet. This is necessarily the case for a local model in the 
case of constant specific vorticity, because of symmetry ar- 
guments. Effects originating a distance ~ H from the planet 
such as the production of the Lindblad wakes can destroy 
this symmetry by means of a back pressure effect. This in- 
creases the gradient of w near the planet, which, if strong 
enough, can give rise to three possible stagnation points. 
For the case shown in Fig. [S] the back pressure is such that 
one stagnation point appears, but shifted in azimuth by ap- 
proximately one softening length. In Fig. [5] we show the 
j/-derivative of w + $l at the radial location of the stagna- 
tion point. 

When no cut off is used, the back pressure gives rise 
to a gradient of w at y = 0, pushing the stagnation point 
away from the planet. The actual configuration of the stag- 
nation points depends on details in the flow, fo r example the 
backg round surface density gradient. Indeed, iMasset et al.l 
(2006) reported three stagnation points for the case with 
constant background surface density, a configuration that 
we find as well for the same surface density gradient. Since 
both the back pressure and the potential are proportional to 
q in the linear regime, the location of the stagnation points 
do es not depend on the planet mass, which was also reported 
bv lMasset etafl i|2006r i. 

When b <C h, the Lindblad torques will be largely in- 
dependent of the softening parameter, since they are gen- 
erated at distances much larger than br p in this case. We 
then expect any related back pressure to be small, with any 
stagnation points located close to the planet. We find this 
indeed to be the case, and we see that for b <C h the agree- 
ment of x s as measured from the simulations with that found 
using the model of Section [3] which does not take account 
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Figure 9. Derivative of the sum of the enthalpy and the potential 
at the orbital radius of the stagnation point for q/h 3 = 0.0252 
and b/h = 0.6, with and without the cut off procedure. 
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Figure 10. Half-width of the horseshoe region vs. softening pa- 
rameter. Two planet masses are considered, the lowest one with 
and without the cut-off procedure. For the more massive planet 
x B /2 is shown, to remove the y/q dependence. Error bars indicate 
an error of 10%. 

of phenomena related to epicyclic motions such as Lindblad 
torques, is much better. 

5.3 Width of the horseshoe region 

We have measured the half- width of the horseshoe region x s 
for various values of b and q though a streamline analysis. 
The results are shown in Fig. 1101 with 10% error bars in- 
dicating our error estimate, for the cases with and without 
applying the cut off procedure. As indicated in our above dis- 



cussion, horseshoe widths are always larger when the cut off 
procedure is applied. The deviation varying from about 10% 
at very small softening to 50% when b ~ h. The results pre- 
dicted by the model in Section [4] with the two-dimensional 
Green's function are indicated by the solid curve. These are 
in good agreement with the other results for small soften- 
ing but give values for x s ~ 20% larger than those found 
from the simulations without a cut off procedure. It appears 
that these results fall below the others, because of the effect 
discussed above that we described as being due to a back 
pressure related to phenomena such as the wake produced at 
a distance ~ 2H/3 from the planet. As expected, for smaller 
values of b this effect is reduced. 

In Fig. [10] we also show results for a planet with a mass 
that is four times larger (corresponding to 4 orbiting a 
Solar mass star, embedded in a disc with h = 0.05). For this 
case, the measured value of :r s was divided by 2 to remove 
the y/q scaling. If x B oc y/q in this mass range, the black 
circles and diamonds should fall on top of each other in Fig. 
1101 It is clear that for all values of the softening parame- 
ter that we cons ider, the hors e shoe width scales as y/q in 
this mass range. iMasset et al.l (|2006r ) speculated that this 
scaling would brake down for softening parameters smaller 
than the Bondi radius. We find this not to be the case, since 
our smallest softening parameter b < q/h 2 for both planet 
masses placing it in that regime. 

Although the perturbed surface density at the location 
of the planet can be quite large for small softening, this per- 
turbation is almost in hydrostatic equilibrium and so does 
not play a major role in the flow, and does not cause a de- 
parture from a; s oc y/q. Below, we will argue that instead, 
while the value of b is important, this departure is governed 
by the ratio q/h? . 

Additional models with different values of H, while 
keeping b/h fixed, confirm that x s oc 1/yh (see equation 
H39|> ) in the same mass range, as long as i s C fl. 

5.4 Horseshoe drag 

IWardl { 199l[ ) found an expression for the corotation torque 
produced by material in the coorbital region in the form 

r Clhs =|Q-a)4Sr^ ; (44) 

where the surface density E oc r~ a , a nd all quantit i es are 
evaluated at r — r p . On the other hand lTanaka et all (|2002l ) 
found an expression for the corotation torque derived from 
linear perturbation analysis in the form 

r c ,iin = 1.36 Q - a\ j£Er*n£. (45) 

Here we shall make use of these expressions, referring the 
reader to a companion paper for additional discussion. 

We first note that although they are both torques, 
these expressions were derived for systems with differing flow 
topology and so should not be expected to be the same. The 
horseshoe drag applies to a system with the flow topology 
of our solutions for the coorbital region which has separatri- 
ces. On the other hand the linear corotation torque is derived 
from the linear perturbation theory of circular orbits. 

The torque in both expressions is proportional to 
the vortensity gradient which cancels out when they are 
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Figure 11. Half-width of the horseshoe region vs. mass ratio q, 
for two different values of b/h. The dashed lines indicate Xslj—g, 
given by equation 1 140 l l. and 0.6 times this value to account for 
the effects of softening and back pressure. The horizontal dotted 
line indicates x B = 2H/3. The tilted dotted line shows equation 
11471 . expected to be valid at high masses. 



equated. In this paper we considered only constant vorten- 
sity for which there is no torque. But we can consider the 
case of a small and smooth vortensity gradient, and follow- 
ing the discussion of section[3]we can argue, as has also been 
confirmed in simulations, that the horseshoe width should 
be c lose to that found as suming constant vortensity. 

iMasset et ail (|2006h assumed that torques obtained 
from simulation results with b/h = 0.3 could be used to de- 
term ine values of x s by e quating them to torques obtained 
from iTanaka et all (|2002h even though the latter were cal- 
culated for 6 = 0. The determined values of a; s agreed with 
those directly measured from the simulations, and we re- 
mark th at our measurement s of x s for b/h = 0.3 agree with 
those of IMasset et all i|200fj) . 

However, from Fig. 1101 we see that x s increases by a 
factor of 1.67 as b decreases from 0.3h to zero. Therefore, 
the horseshoe drag torque, being proportional to 
tually nearly an order of magnitude larger than the linear 
corotation torque for 6 = and so should not be equated to 
it. It thus turns out that the combined effects of finite soft- 
ening and the back pressure phenomenon discussed above 
reduced the estimated horseshoe drag by almost an order 
of magnitude compared t o the value that sho uld have been 
adopted to compare with ITanaka et al. I <|2002l ). 



5.5 Extension to higher masses 

When q/h 3 is of order unity, the Hill sphere, the Bondi ra- 
dius and the half-width of the horseshoe region are all com- 
parable to H . At approximately this mass, the waves excited 
by the planet start to become non-linear at a distance ~ H 
from it, and gap formation sets in. This reduces the strength 
of the Lindblad torques and any influence of material at a 
distance ~ H, and therefore the back pressure effect de- 



scribed above. At the same time, the width of the horseshoe 
region, being proportional to \J q/h, also becomes compara- 
ble to H, extending into the Lindblad resonance region. This 
also reduces the back pressure effect, making the horseshoe 
width larger. This was checked by running an additional 
model with q/h 3 = 0.4032 (8 embedded in a disc with 
h = 0.05 around a Solar mass star), b/h = 1 and b/h — 0.2. 
For b/h = 1, we find a horseshoe width of x s = 0.54//, which 
corresponds to an exact scaling x s oc y/q according to Fig. 
1101 This is to be expected, since x s < 2Z//3, and the horse- 
shoe region does not extend into the wave excitation zone. 
For b/h = 0.2,, we find a; s — 0.9//, while a scaling with ^/q 
would imply x s = 0.75//, according to Fig. [10] Therefore, 
the horseshoe width is larger than predicted by the sim- 
ple m odel, in agreement with the findings of IMasset et al.l 
l|2006h . We stress that this behaviour should not be seen 
as an onset of nonlinearity, since there is no horseshoe re- 
gion in linear theory, but rather as a reduction of the back 
pressure effect of material at a distance ~ 2Z//3 from the 
planet where the main Lindblad torques are produced on 
the horseshoe region. 

The exact mass at which x s = 2Z//3 of course de- 
pends on softening. If we write x s = C(b)x s (b — 0), with 
C(b) < 1 and for an approximate estimate, set x B (b — 
0)/r p = 1.68\/ q/h (see equation (|40|l ). we get: 



q < Qtrana 



0.157/t 3 

c(by 



(46) 



as a condition for x s oc ^fq. For b/h = 0.3, as used by 
IMasset et aD (|2006T ). we have C w 0.6 (see Fig. HDJ, and 
therefore gtrans = 5.5 ■ 10~ 5 for h = 0.05, exactly the 
mass at which the departure from linearity as reported by 
IMasset et al. | i|200rj) be gins. Note that for the same param- 
eters, the ratio of the Bondi radius to the softening param- 
eter is approximately unity. However, for C = 1, we have 
qtrans/h 3 = 0.157, confirming that for q/h 3 = 0.1008 and 
q/h 3 = 0.0252, as shown in Fig. 1101 x B oc ^fq for all soften- 
ing parameters, independent of the ratio of the Bondi radius 
to the softening. On the other hand, for b^> h, using equa- 
tion (I41|) . it can be shown that qtrnns/h 2 oc b. In this case, 
the critical value of b for which q — gtrans is proportional to 
the Bondi radius. However, such a large softening is incom- 
patible with the idea of vertical averaging, from which we 
expect b to be of the order of h. 

These ideas are further illustrated in Fig. 1111 where we 
show the measured half-width of the horseshoe region ver- 
sus planetary mass. For q <C h 3 , we expect x B oc y/q, with 
a coefficient that depends on the softening. For small soft- 
ening, we see that the results agree very well with equation 
(|40|) , while for large softening we find values that are a factor 
of 0.6 smaller due to the combined effect of non-zero soft- 
ening and back pressure (see above). At higher masses, we 
expect x B to be proportional to the radius of the Hill sphere 
(Pepliriski 2008): 



2.47 



(I) 



1/3 



(47) 



indicated by the tilted dotted line in Fig. 1111 For q ^> h we 
find this indeed to be the case, independent of the value of 
b/h. The results for b/h = 0.4 are in very good agreement 
with IMasset et all ([2006) , who used b/h = 0.3, for all values 
of q. In between the two regimes of x s oc ^fq and x s oc q 1 ^ 3 
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the width of the horseshoe region rises faster than Jq, which 
is due to a reduction in the strength of the back pressure, as 
argued above. For b/h = 0.05, no such behaviour is found, 
since the effects of the back pressure are small for all q. 

We therefore conclude that for q > gtrans, the back pres- 
sure effect of the Lindblad torques is reduced, which leads to 
an increase in the width of the horseshoe region towards the 
value obtained from equation (|39[1 . T his is consistent wit h 
the streamline analysis presented in iMasset et al. I (|2006h . 
where it is shown that in this transition regime only one 
stagnation point survives and moves towards the location of 
the planet. We have confirmed this behaviour in our simu- 
lations. This increase in x s can have a major impact on the 
torque acting on the planet (see equation (|44[l ). 



6 DISCUSSION AND CONCLUSIONS 

In this paper we presented a simple model of the coorbital 
region around a low mass planet. Using this we derived the 
horseshoe width as a function of planet mass and gravita- 
tional softening parameter. In the limit of zero softening we 
found that 



(48) 



This result agreed with high resolution numerical simu- 
lations to within several percent. However for softening 
lengths br p ~ H, the discrepancy was larger, with the simu- 
lations indicating a horseshoe width about 22% smaller. By 
considering simulations for which a cut off procedure was 
used to remove the effects of the protoplanet potential pro- 
duced at and beyond a radial separation of 2/3H from it, it 
was found that phenomena at that separation could signifi- 
cantly affect the horseshoe width, even when that was much 
narrower, distorting the streamlines and reducing the width 
through the action of an additional back pressure that is 
more effective for larger softening. This may artificially re- 
duce the horseshoe drag in such cases. 

We also used our results to show that the horseshoe 
drag, exerted by material executing horseshoe turns is about 
an order of magnitude larger than the linear corotation 
torque in the zero softening limit. A more complete compar- 
ison between linear corotation torques and horseshoe drag 
for finite b requires additional linear calculations which are 
presented in detail in an accompanying paper. There we also 
find that the non-linear corotation torque (horseshoe drag) 
is always much larger than the linear corotation torque for 
non zero b. 

We have focused on a two-dimensional description of 
the horseshoe region, with a softening parameter b in the 
planet potential which may approxi mately account for t hree- 
dimensional effects. As reported in lMasset et al.l l|2006l ). the 
horseshoe drag torque appears to be stronger in fully three- 
dimensional simulations compared to two-dimensional runs 
that include softening. Clearly, a three-dimensional model 
of the horseshoe region is desirable. This will be the subject 
of a future investigation. 

Another useful extension of the present discussion 
would be the inclusion of non-barotropic effects. We remark 
that the model presented in this paper is valid for discs that 
have a constant specific vorticity and entropy, the latter con- 
dition leading to a barotropic equation of state. Introducing 



a radial vortensity gradient breaks the up-down symmetry 
in Figs. 1-3, but simulations show that this effect is barely 
detectable. Thus we may expect that the prediction of the 
horseshoe width obtained from our simple model may work 
reasonably for non-barotropic discs with a radial entropy 
gradient. However, we do expect some difference in x s be- 
tween isothermal and non-isothermal discs. It is easy to see 
that obtained from equation (136[1 . is proportional to 

c^ 1 ^ 2 . For equal temperatures, the sound speed in an adia- 
batic disc is a factor s ffi larger than the isothermal sound 
speed, where 7 is the adiabatic exponent. This makes the 
horseshoe region a factor 7 1 / 4 smaller in adiabatic simu- 
lations. Although the difference lies within our 10% error 
bars, we have no ticed it when comparing our prese nt re- 
sults with those in lPaardekooper fc Papaloizou! l|2008f ). Note 
that this makes the adiabatic horseshoe drag, being propor- 
tional to a;*, a factor 7 smaller than the isothermal horseshoe 
drag. Since also the wave torque scales as -y~ , the relative 
strength of the Lindblad torques and horseshoe drag remains 
the same for adiabatic discs. 

The shape of the horseshoe region changes when a 
global radial mass flow is introduced with respect to the 
planet. This mass flow can be due to viscous accretion, but 
also due to radial movement of the planet when allowing the 
orbit of the planet to change. When the time scale of the ra- 
dial flow with respect to the planet to cross the horseshoe 
region is smaller than the libration time scale, an asymmetry 
between the sides of the horseshoe region l eading and trail- 
ing the planet develops (|Artvmowic 3 l2004l). This is impor- 



tant for studying Type III migration ( Masset fc Papaloizou! 
2003). In this paper, we have kept the planet on a fixed orbit 
in an inviscid disc, and therefore such effects did not occur. 
More work is necessary to study the importance of including 
the effect of planetary migration on the disc response and 
torques for low-mass planets. 
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